R Notebook

• Load necessary libraries

library(ggplot2)
library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(plotly)

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

Attaching package: ‘plotly’

The following object is masked from ‘package:ggplot2’:

    last_plot

The following object is masked from ‘package:stats’:

    filter

The following object is masked from ‘package:graphics’:

    layout

library(scales)

• Load the dataset

df <- read.csv("Forest/forest.csv")
head(df)

• Check for any null values

any(is.na(df))

[1] FALSE

Data Cleaning

• Combine all the soil types into a single column

# Since there are 40 of the columns dedicated only to soil type, it was necessary to reduce those columns into one
df <- df %>%
  mutate(Combined_Soil_Type = case_when(
    Soil_Type1 == 1 ~ "1",
    Soil_Type2 == 1 ~ "2",
    Soil_Type3 == 1 ~ "3",
    Soil_Type4 == 1 ~ "4",
    Soil_Type5 == 1 ~ "5",
    Soil_Type6 == 1 ~ "6",
    Soil_Type7 == 1 ~ "7",
    Soil_Type8 == 1 ~ "8",
    Soil_Type9 == 1 ~ "9",
    Soil_Type10 == 1 ~ "10",
    Soil_Type11 == 1 ~ "11",
    Soil_Type12 == 1 ~ "12",
    Soil_Type13 == 1 ~ "13",
    Soil_Type14 == 1 ~ "14",
    Soil_Type15 == 1 ~ "15",
    Soil_Type16 == 1 ~ "16",
    Soil_Type17 == 1 ~ "17",
    Soil_Type18 == 1 ~ "18",
    Soil_Type19 == 1 ~ "19",
    Soil_Type20 == 1 ~ "20",
    Soil_Type21 == 1 ~ "21",
    Soil_Type22 == 1 ~ "22",
    Soil_Type23 == 1 ~ "23",
    Soil_Type24 == 1 ~ "24",
    Soil_Type25 == 1 ~ "25",
    Soil_Type26 == 1 ~ "26",
    Soil_Type27 == 1 ~ "27",
    Soil_Type28 == 1 ~ "28",
    Soil_Type29 == 1 ~ "29",
    Soil_Type30 == 1 ~ "30",
    Soil_Type31 == 1 ~ "31",
    Soil_Type32 == 1 ~ "32",
    Soil_Type33 == 1 ~ "33",
    Soil_Type34 == 1 ~ "34",
    Soil_Type35 == 1 ~ "35",
    Soil_Type36 == 1 ~ "36",
    Soil_Type37 == 1 ~ "37",
    Soil_Type38 == 1 ~ "38",
    Soil_Type39 == 1 ~ "39",
    Soil_Type40 == 1 ~ "40"
  ))

• Doing the same for wilderness area

df <- df %>%
  mutate(Combined_Wilderness_Area = case_when(
    Wilderness_Area1 == 1 ~ "Rawah Area",
    Wilderness_Area2 == 1 ~ "Neota Area",
    Wilderness_Area3 == 1 ~ "Comanche Peak Area",
    Wilderness_Area4 == 1 ~ "Cache la Poudre Area"
  ))

• Summing up all hillshade into one

df <- df %>%
  mutate(Total_Hillshade = Hillshade_9am + Hillshade_Noon + Hillshade_3pm)

• Change degrees into cardinal directions

# Function to convert Aspect into cardinal and intermediate directions
get_cardinal_direction_full <- function(aspect) {
  if (!is.na(aspect) && aspect >= 0 && aspect < 360) {
    if (aspect >= 337.5 || aspect < 22.5) {
      return("N")
    } else if (aspect >= 22.5 && aspect < 67.5) {
      return("NE")
    } else if (aspect >= 67.5 && aspect < 112.5) {
      return("E")
    } else if (aspect >= 112.5 && aspect < 157.5) {
      return("SE")
    } else if (aspect >= 157.5 && aspect < 202.5) {
      return("S")
    } else if (aspect >= 202.5 && aspect < 247.5) {
      return("SW")
    } else if (aspect >= 247.5 && aspect < 292.5) {
      return("W")
    } else if (aspect >= 292.5 && aspect < 337.5) {
      return("NW")
    }
  }
  return("N") 
}

# Apply the function to the Aspect column directly
df$Aspect <- sapply(df$Aspect, get_cardinal_direction_full)

# View the first few rows to confirm
head(df$Aspect)

[1] "NE" "NE" "SE" "SE" "NE" "SE"

• Set the values for cover type

cover_type_mapping <- c(
  `1` = "Spruce/Fir",
  `2` = "Lodgepole Pine",
  `3` = "Ponderosa Pine",
  `4` = "Cottonwood/Willow",
  `5` = "Aspen",
  `6` = "Douglas-fir",
  `7` = "Krummholz"
)

# Replace the Cover_Type values in the dataframe
df$Cover_Type <- factor(cover_type_mapping[as.character(df$Cover_Type)], levels = cover_type_mapping)

# View the first few rows to confirm the changes
head(df$Cover_Type)

[1] Aspen          Aspen          Lodgepole Pine Lodgepole Pine Aspen          Lodgepole Pine
Levels: Spruce/Fir Lodgepole Pine Ponderosa Pine Cottonwood/Willow Aspen Douglas-fir Krummholz

• Since the Slope ranges from 0 to 50, we can make it into continuous value so that while plotting it makes it easier and not very messy

# Modify the Slope column to group values in bins of 5
df$Slope <- cut(df$Slope, 
                breaks = seq(0, max(df$Slope, na.rm = TRUE), by = 5), 
                include.lowest = TRUE, 
                right = FALSE, 
                labels = paste(seq(0, max(df$Slope, na.rm = TRUE) - 5, by = 5), 
                               seq(5, max(df$Slope, na.rm = TRUE), by = 5), sep = "-"))

• Check distribution for Elevation

a <- ggplot(df, aes(Elevation)) +
  geom_histogram(fill="steelblue")+
  labs(title= "Distribution of Elevvation")
a

a+facet_wrap(Combined_Wilderness_Area~.) +
   labs(title = "Wilderness Areas across different elevation levels")

ggplot(df, aes(Elevation)) +
  geom_histogram(fill="steelblue")+
  labs(title= "Distribution of Elevation") +
  facet_wrap(Aspect~.)

ggplot(df, aes(x = Elevation)) +
  geom_density(fill = "blue", alpha = 0.5) +
  labs(title = "Density Plot of Elevation", x = "Elevation", y = "Density") +
  theme_minimal()

• Since our data set is too large for proper visualization, we will take a small sample of the data that resembles the trend of the original data.

set.seed(123)  # Setting a seed for reproducibility
df_sample <- df %>% sample_n(3000)

ggplot(df_sample, aes(x = Elevation)) +
  geom_density(fill = "blue", alpha = 0.5) +
  labs(title = "Density Plot of Elevation", x = "Elevation", y = "Density") +
  theme_minimal()

• Check the distribution of Slope

ggplot(df, aes(Slope)) +
  geom_bar()

ggplot(df, aes(x = Slope)) +
  geom_density(fill = "blue", alpha = 0.5) +
  labs(title = "Density Plot of Slope", x = "Slope", y = "Density") +
  theme_minimal()

ggplot(df_sample, aes(x = Slope)) +
  geom_density(fill = "blue", alpha = 0.5) +
  labs(title = "Density Plot of Slope", x = "Slope", y = "Density") +
  theme_minimal()

Since the density of slope and elevation looks similar to that of the real dataset, we will conclude that other feature will follow the same trend.

• Distribution of Cover Type throughout the data set

a<-ggplot(df, aes(Cover_Type)) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_y_continuous(labels = scales::comma)+
  geom_bar(fill="steelblue")+
  geom_text(stat = 'count', aes(label = ..count..), 
            vjust = -0.5, color = "black") +
  labs(title = "Distribution of Cover Type",
       x="Cover Type",
       y = "Count")
a

• Facet the cover distribution according to Wilderness Area

a + facet_grid(Combined_Wilderness_Area~., scales="free")

a+facet_wrap(Aspect~.)

• Check how much sunlight does each direction receive at 9am

ggplot(df, aes(Hillshade_9am))+
  geom_histogram(fill="steelblue") +
facet_wrap(Aspect~.) +
  labs(title = "Hillshade at 9am")

East and North received the most sunlight at 9am.

• Check how much sunlight does each direction receive during noon

ggplot(df, aes(Hillshade_Noon))+
  geom_histogram(fill="steelblue") +
  facet_wrap(Aspect~.) +
  labs(title = "Hillshade during Noon")

During the noon time, only south west and west received the least amount of sunlight.

• Check how much sunlight does each direction receive at 3pm

ggplot(df, aes(Hillshade_3pm))+
  geom_histogram(fill="steelblue") +
  facet_wrap(Aspect~.) +
  labs(title = "Hillshade at 3pm")

From all the plots above, it is clear that west and southwest received the least amount of sunlight.

• Remove the unnecessary columns Since our data set already had a massive set of columns and the columns ‘Soil_Type’ and ‘Wilderness_Area’ have been categorized into the column ‘Combined_Soil_Type’ and ‘Combined_Wilderness_Area’ they are unnecessary

df <- df %>%
  select(-starts_with("Soil_Type"), -starts_with("Wilderness_Area"), -starts_with("Hillshade")) #Remove any column that starts with the word Soil_Type, Wilderness_Area and Hillshade

• Checking to see if the changes have been applied

df

a+facet_wrap(.~Combined_Wilderness_Area)

• Plot density graph for elevation vs cover type

ggplot(df, aes(x = Elevation, fill = factor(Cover_Type))) +
  geom_density(alpha = 0.5) +
  labs(title = "Density of Elevation by Cover Type", x = "Elevation", y = "Density")

We can see that Cottonwood grows at lower altitude and Krummholz grows at the higher altitude.

ggplot(df_sample, aes(x=Cover_Type, y=Elevation))+geom_boxplot()

• Check distribution for soil type

p<-ggplot(df, aes(x = Combined_Soil_Type)) +
  geom_bar(fill = "steelblue") +
  labs(title = "Distribution of Combined Soil Types", x = "Soil Type", y = "Count") +
  theme(axis.text.x = element_text(angle = 0, hjust = 1))
p

From the bar plot, it can be seen that the most type of prevalent soil is of type 29 which is almost 1200000. On the other hand, Most of the soil types are almost non-existent like soil_Type 13,14,37 and many more which indicates their presence in somewhat small amount compared to others.

• Facet soil type distribution by wilderness area

p+facet_wrap(Combined_Wilderness_Area~.) +
  labs(title = "Distribution of Soil Type across Wilderness Areas")

Rawah Area receives has the most amount of soils while Comanache Area has the most variety of soils.

•

ggplot(df, aes(x = Cover_Type, y = Total_Hillshade, fill = Cover_Type)) +
  geom_boxplot() +
  labs(title = "Total Hillshade (Sunlight) Distribution by Cover Type",
       x = "Cover Type",
       y = "Total Hillshade (Sunlight)") + scale_y_continuous(labels = scales::comma) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) + facet_grid(Combined_Wilderness_Area~.)

From the graph above, it seems that every type of plant can grow when receiving more sunlight but Lodgepole Pine can grow on areas with low sunlight of upto only 200 units.

• Create a bar for soil type and stack it with cover type

ggplot(df, aes(x = Combined_Soil_Type, fill = Cover_Type)) +
  geom_bar(position = "stack") +
  labs(title = "Cover Type and Soil Type", x = "Soil Type", y = "Count") +
  theme(axis.text.x = element_text(angle = 0, hjust = 1))

From the bar plot we can see that soil 29 is the most suitable for spruce/fir tree and Lodgepole Pine tree which are also the most present type of trees on other soils compared to other plants indicating their capacity to grow on any type of situation.

• Facet it with Wilderness Area

ggplot(df, aes(x = Combined_Soil_Type, fill = Cover_Type)) +
  geom_bar(position = "stack") +
  labs(title = "Cover Type and Soil Type", x = "Soil Type", y = "Count") +
  facet_wrap(Combined_Wilderness_Area~.)+
  theme(axis.text.x = element_text(angle = 0, hjust = 1))

• Check to see if the relation between soil type and hydrology

b<-ggplot(df, aes(x=Horizontal_Distance_To_Hydrology, y=Vertical_Distance_To_Hydrology)) + geom_smooth() + facet_wrap(~Combined_Soil_Type, scales="free") +
  labs(title = "Type of Soil available near water bodies")

From the plots it is hard to determine relation between the two as every soil type seem to be available at almost every horizontal distance away from a body of water. Some like type 1, 16, 18, 2, 21, 8, 9, 5, 6, 7 is found in the range of 600m away from the body of water, the other types are available 1000 meters away. As for the vertical distance from hydrology, most are found at the range of 100-150 meters above and some like 8 and 25 are found below the body of water. Soil type 4, 24, 27, 28, 31, 32, 38, 39 and 40 could be found above 200m above any body of water but soil type like 7, 8, 9 and 25 are found as close as 40 meter and even potrayed in negative indicating them being submerged in the body of water

• Add new feature ‘Cover_Type’ as colour

ggplot(df_sample, aes(x=Horizontal_Distance_To_Hydrology, y=Vertical_Distance_To_Hydrology, colour = Cover_Type)) + geom_point() + facet_wrap(~Combined_Soil_Type, scales="free") +
  labs(title = "Cover type and sources of water")

• Point graph with Elevation vs Soil Type

ggplot(df_sample, aes(y=Elevation, x=Combined_Soil_Type, colour = Cover_Type))+
  geom_point()+
  labs(title = "Elevation vs Soil Type")

As the type of soil increases so does the elevation.

•

# Create a boxplot
ggplot(df_sample, aes(x = Elevation, y = Total_Hillshade, colour = Cover_Type)) +
  geom_point(alpha = 0.5) +
  labs(title = "Sunlight Exposure (Hillshade) vs Elevation",
       x = "Elevation",
       y = "Total Hillshade (Sunlight)") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  facet_grid(Combined_Wilderness_Area~.)

• Check the

ggplot(df_sample, aes(x = Elevation, y = Total_Hillshade)) +
  geom_point(alpha = 0.5) +
  labs(title = "Sunlight Exposure (Hillshade) vs Elevation",
       x = "Elevation",
       y = "Total Hillshade (Sunlight)") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

•

ggplot(df, aes(y=Total_Hillshade,x =Cover_Type, fill = Cover_Type)) +
  geom_boxplot() +
  theme(axis.text.x = element_text(angle = 45))

• Plot a histogram forhillshade and facet it with wilderness area and cover type

ggplot(df, aes(x = Total_Hillshade)) +
  geom_histogram(position = "dodge", binwidth = 10, fill="steelblue") +
  labs(title = "Distribution of Cover_Type", x = "Hill Shade(Sunlight)", y = "Count") +
  theme(axis.text.x = element_text( hjust = 1)) + 
  facet_grid(Cover_Type~Combined_Wilderness_Area, scales = "free")

From the plot, we can see that some wilderness area is unsuitable for some type of forest as they do not grow there at all. Eg, Spruce and …… do not grow on wilderness area 4 at all but ….wood only grows on wilderness area 4. If we were to look at the plot alone it would seem like Wilderness_Area4 hosts the most forest but if we look at the count then we can see that Wilderness_Area1 hosts the most forest. If we were to determine the conditions of the Wilderness Area then we can determine the type of forest that grows on such terrains.

q <- ggplot(df, aes(x = Combined_Wilderness_Area, fill = Cover_Type)) +
  geom_bar(position = "stack") +
  labs(title = "Distribution of Cover Type according to Wilderness Area",
       x = "Wilderness Area",
       y = "Count") +
  geom_text(stat = 'count', aes(label = ..count..), position = position_stack(vjust = 0.5), size = 3, color = "white") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
q

ggplotly(q)

From the bar graph, we can see that the most prevalent type of forest is Lodgepole Pine followed by Spruce/Fir which are at an noticeable amount followed by Ponderosa Pine and so on. We can also see that the Rawah Area hosts the most number of plants while the Comanche Peak Area hosts the most type of plants.

• Create a point plot for elevation vs horizontal distance to roadways

a<-ggplot(df_sample, aes(y=Elevation, x=Horizontal_Distance_To_Roadways, size = Slope, colour = Cover_Type)) + geom_point(alpha = 0.5) + labs(title = "Distribution of Cover types near roadways")
a

• Facet the above plot by wilderness area

a+facet_wrap(Combined_Wilderness_Area~., scales = "free")

ggplotly(a)

Warning: Using size for a discrete variable is not advised.

• Create a histogram for Horizontal_Distance_To_Fire_Points and facet it with Cover_Type and Combined_Wilderness_Area

ggplot(df_sample, aes(Horizontal_Distance_To_Fire_Points)) + geom_histogram() +
  facet_grid(Combined_Wilderness_Area~Cover_Type) +
  labs(title="Areas prone to forest fire")

From the graph above, it seems that Rawah Area and Comanache Peak Area are prone to wild fires. Specially species of plant like Spruce and Lodgepole Pine are planted near those areas rather than any other species.

ggplot(df_sample, aes(x=Horizontal_Distance_To_Fire_Points,y=Horizontal_Distance_To_Hydrology, colour = Cover_Type)) + geom_point() +
  facet_wrap(.~Combined_Wilderness_Area)

ggplot(df_sample, aes(x=Horizontal_Distance_To_Roadways, y=Horizontal_Distance_To_Fire_Points, colour = Cover_Type)) +
  geom_point() +facet_wrap(Combined_Wilderness_Area~.)

From the graph it can be seen that in Cache la Poudre Area is closer to Roadways and Fireplaces but the plants are at a scare amount nut in Neota Area, Spruce/Fir are closest to roadways and fire points but also at an scarce amount. On the other hand, in Comanche Peak Area Spruce/Fir, Lodgepole pine are present at a high amount and are closer to roadways and fire point but Aspen and Douglas fir like plants are planted further awat from those points. The same could be said for the Rawah Area except for this area Spruce/Fir and Aspen are planted close to fire points.

ggplot(df_sample, aes(x=Horizontal_Distance_To_Hydrology, y=Horizontal_Distance_To_Fire_Points)) +
  geom_point()

https://roh4n.shinyapps.io/Rohan_Kalu_23189640/